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Abstract. In this paper we propose a framework for the analysis of the 
cardiac and respiratory motion of the myocardium based on stereo video 
images taken from a surgical robot. Our aim is to use the analysis results 
as the basis for a myocardial surface reconstruction algorithm. These 
results will provide important prior information that will be needed to 
obtain accurate surface point matching in homogeneous regions and in 
areas of high deformation. 

We show how to fit a parametric cardio-respiratory model to myocardial 
surface features, and we present methods for validating the model. 


1 Introduction 


There have been a number of recent technological advances in cardiac surgery 
aimed at making the procedures less invasive. These include minimally-invasive 
direct coronary artery bypass (MIDCAB) and totally endoscopic coronary artery 
bypass (TECAB). In TECAB, the entire procedure is performed robotically 
through small incisions on the beating heart, without the need for a cardiopul- 
monary bypass. 

However, TECAB procedures still suffer from a relatively high rate of con- 
version to either MIDCAB or conventional open-chest surgery [1-3]. Reasons 
for this include difficulty in locating an artery or vessel misidentification. Aug- 
mented reality (AR) guidance has the potential to reduce the conversion rate 
by overlaying an accurately registered graphical representation of a preoperative 
imaging model onto the surgeon’s view of the patient. 

The da Vinci? robot provides an ideal device for such guidance due to 
the high quality 3D visualisation it provides. Previous researchers have demon- 
strated AR on phantom or animal data [4-6] or visualisation of intraoperative 
ultrasound [7]. Our aim is to augment the da Vinci?™ robot’s display with a 4D 
model of the beating heart preoperatively constructed from clinical data. 

In order to register a preoperative 4D model to the robot’s intraoperative, 
stereo video images, it is necessary to reconstruct the 4D myocardial surface from 
the stereo video images. In [8], Stoyanov presents a method for sparse myocardial 
surface reconstruction based on tracking features in 3D, and in [9] he describes 
a method for dense reconstruction based on piecewise bilinear maps. Neither 
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of these methods take into account prior information about myocardial motion 
however (apart from assuming smoothness). Such information could potentially 
be used to improve point matching in areas of high deformation, and to infer 
more accurate interpolation methods for reconstructing the large homogeneous 
regions in which point matching fails. 

The value of exploiting the periodicity of myocardial motion is shown by Ort- 
maier in [10], where the trajectory of an occluded feature is accurately predicted 
based on its past motion and the motion of other visible features, with which its 
motion is correlated. Accurate predictions such as this can be incorporated into 
a motion model to provide strong prior information for feature tracking. 

In this paper we will focus on measuring the myocardial motion and decom- 
posing it into cardiac and respiratory components. Knowing how much cardiac 
and respiratory motion is likely to occur between two consecutive video frames 
provides us with prior information that allows us to put empirically-supported 
bounds on the state space that a feature tracking algorithm must search. In or- 
der to achieve accurate sparse 4D reconstruction (which will be the initial stage 
of our dense reconstruction), a tracking algorithm must make effective use of 
such prior information. 


2 Method 


Our experiment is based on the trajectories of myocardial features manually 
tracked in stereo intraoperative videos of a patient. The trajectories were re- 
constructed in 3D after calibrating the cameras. Before we could measure the 
cardiac and respiratory components of the motion of the features, we had to 
decompose the trajectories into these components. We based our decomposition 
on the following assumptions: 


— The motion of each feature is the sum of cardiac motion, respiratory motion 
and some amount of noise (noise being a result of tracking errors and random 
fluctuations in the cardio-respiratory motion). 

— The cardiac and respiratory motions are periodic up to the effects of noise’. 


2.1 Motion Decomposition 


In [11] a Cardiac Respiratory Parametric Model (CRPM) is used to describe the 
motion of the coronary arteries extracted from angiographic images. We have 
adopted this model to characterize the motion patterns of myocardial surface 
features. The CRPM models the motion of a myocardial feature as the sum of 
two periodic B-splines with uniform knot vectors, parameterised by the cardiac 
phase and the respiratory phase. 


' The patients are under general anaesthetic during the operation, so variations in 
cardiac motion should be small. It is also reasonable to assume that variations in 
the respiratory motion are small, as the patient’s breathing is mechanically controlled 
by a respirator. 


So, suppose the trajectory of feature f in an arbitrary dimension is given 
by: wf = (x f0(X0,P0),---,@4,7-1(X7-1, PT-1))", where x; € [0,1) and p; € 
(0,1) are measures of the cardiac and respiratory phases at time t. The CRPM 
approximates each x ¢,4(x+, pr) as Ty4(V, Xt, Pr), defined as follows: 
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where: p, /p, are the degrees of the cardiac/respiratory B-splines; n,/n, are the 
numbers of cardiac/respiratory B-spline control points; v = (vo,.--,Un,—p,y—1; 
Unys ee fe ee is a vector of the first ny — py cardiac B-spline control 
points followed by the first n, — pp respiratory control points (the last p,/p, 
control points of the cardiac/respiratory B-splines are a repetition of the first 
py/Pp. This makes the B-splines periodic); b©, and bf, are the degree d basis 
functions at control point 7 for the knot vectors associated with the cardiac and 
respiratory B-splines respectively. 

We want to find the control points that minimise the squared error between 
the estimate and the true feature trajectory: 
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The linear least squares solution to (2) is given by the product of the pseu- 
doinverse of a matrix B of the blending coefficients and «+z. We used a Tikhonov 
regularisation matrix S' of second-order central difference coefficients to minimise 
the curvature of the respiratory component to a degree controlled by a. 


v=(B'B+a07S'S) Ba; , (3) 


The reader is referred to [11] for a more detailed description of the terms of this 
solution. 


2.2 Parameter Estimation by k-Fold Cross Validation 


As in [11], we used cubic B-splines (p, = pp, = 3). We used k-fold cross validation 
to find values for ny, n, and a in a finite search space that would: (7) minimise the 
model’s error over test data that was omitted from the model fitting procedure; 
and (ii) minimise the Cardiac Cycle Registration (CCR) error (see next section) 
over the test data. 

Given values for ny, N,) and a, we calculated the average prediction error 
e(ny, Np, a) as follows: 


— Partition the trajectory af into k equal parts, omit part i and calculate (3) 
with the remaining k — 1 parts. 


— Evaluate the RMS error \/£;[e?] between the model and the omitted part, 
and the CCR error €€°R over the omitted part. 
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— e(ny, Np, a ay the sum of the RMS model fitting errors 
and the CCR errors, averaged over the k parts. 


Large values of k are preferable, as they allow e(n,,n,,a) to be calculated 
over more test data sets. However the value of k chosen should be small enough 
for each partition to contain at least one respiratory cycle. 


2.3. Cardiac Cycle Registration 


No motion decomposition validation was given in [11]. However under the peri- 
odicity assumptions stated above, subtracting the respiratory component of the 
CRPM from the data should leave us with the periodic cardiac motion and a 
small amount of noise. Hence, after subtracting the respiratory component and 
temporally registering the cardiac cycles of the resultant data to each other, we 
should be left with a set of cardiac trajectories that are accurately registered 
both temporally and spatially. The spatial variance in the temporally-registered 
trajectories thus gives us a measure of the accuracy of the CRPM’s respiratory 
motion estimation. 

To make this more concrete, suppose the CRPM is parameterised with control 
point vector v. We remove the model’s respiratory component from the data with 
vp, = Xpi(Xir Pi) — Tf i(v,0,p:). We then draw m samples at corresponding 
phases from each cardiac cycle of «4 (uniformly-spaced in their phases). Let 
these samples form an m x c matrix X, where: X;,; is the i*" sample of the j™ 
cycle and c is the number of cardiac cycles in we. The CCR error €C°® is taken 
as the standard deviation of samples drawn from a given cardiac phase, averaged 


over all sampled phases, i.e.: 
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2.4 Estimating Cardiac and Respiratory Phase 


Before we can make use of the above methods, we need an estimate of the cardiac 
and respiratory phase at each tracked point in the feature trajectory. It would 
be ideal if we could acquire this information from the anaesthetist’s monitor, 
but not all such monitors provide an interface through which this information 
can be streamed. Furthermore, monitors that do support such streaming may 
not be exactly in phase with the intraoperative videos. So as an alternative, we 
propose the following phase estimation method based on the autocorrelation of 
the feature trajectories. 

Let xf be the it” element of « f, with unknown cardiac and respiratory 
phases x; and p;. We estimate the cardiac periods of all cycles relative to the 


The trajectories of 5 features over the first cardiac cycle 
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Fig. 1. (Left) The first frame of the video with the 5 tracked features. (Right) A graph 
of the trajectories of the 5 features over the first cardiac cycle (covering about 1 second 
of video). The markers indicate some corresponding frames. The circular markers are 
at the first frame, and the square markers at the last. 


phase of x fo by finding the temporal offsets A that locally maximise the following 
autocorrelation function: 
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where Zf,a:) and Oz;,a:) are the mean and standard deviation of tf.a, Tf,a41; 
Se 

If Ao and A, are consecutive peaks in this function, we estimate the cardiac 
phase x; of the intermediate point t € [Ao, A1) by linear interpolation: 


ula. Cpr : (6) 


To estimate respiratory phase, we first of all filter the high frequency cardiac 
motion out of «» with a Gaussian convolution kernel of a suitable standard devi- 
ation (we used a standard deviation equal to a third of the mean cardiac period), 
and then we use the same autocorrelation method on the filtered data. Note that 
the precision of the phase estimates obtained by this method is limited by the 
frame rate of the video. To obtain higher-precision estimates, we supersampled 
xy using spline interpolation, before calculating fauto. 


3 Results and Discussion 


We manually tracked 5 features over 1060 frames of a 50Hz video of a single pa- 
tient. For given values of n,, mp) and a, we fitted the CRPM to each dimension 
of each feature individually, and then used 4-fold cross validation (as described 


Table 1. This table summarises the RMS model fitting errors and the mean heights 
of the CCR error bars (see figure 2) in each dimension for the CRPMs that were fitted 
to the five features with parameters n, = 19, n, = 8 and a = 1.42. The percentages in 
the ‘mean CCR error bar heights’ columns give the ratios of the mean heights of the 
red error bars in figure 2 to the mean heights of the green error bars. The error bounds 
in the ‘RMS model fitting errors’ columns are the standard deviations of the absolute 
errors. 


Mean CCR error bar heights RMS model fitting errors 
Feature |x (mm) y (mm) z (mm) z(mm) y(mm)~— z (mm) 
1 4.91 (32%) 4.61 (35%) 11.41 (90%) |2.1941.34 1.9241.15 4.7143.01 
2 3.47 (26%) 2.50 (26%) 5.35 (96%) |1.3840.91 1.05+0.66 1.89+1.12 
3 3.72 (26%) 2.89 (28%) 4.13 (93%) |1.4940.96 1.10+0.66 1.49+0.95 
4 3.46 (25%) 2.48 (28%) 3.68 (83%) |1.4540.97 1.00+0.62 1.40+0.90 
5 4.39 (28%) 3.31 (28%) 9.65 (90%) |1.8841.14 1.38+0.83 3.7842.38 
Mean: |3.99 (27%) 3.16 (30%) 6.84 (90%) |1.7141.10 1.34+0.86 2.97+2.15 


in section 2.2) to calculate an average value for e(ny,n,,a) over all dimen- 
sions and features. We searched for values of e’s parameters that would locally- 
optimise this average with a hill climbing strategy over the domain (n,,n», a) € 
{5,...,20} x {5,...,20} x {0,0.01,0.02,...,2}. The values we found are given 
in the table and figures. 

In table 1, the mean CCR error bar heights and RMS fitting errors in the z 
dimension are almost twice the errors in the x and y dimensions. This is to be 
expected, as the stereo endoscopic images we worked with have a very narrow 
baseline, and so the features would need to have been matched very accurately 
to attain accurate depth estimation. A good stereo matching algorithm with 
subpixel accuracy could potentially improve the depth estimation. 

The table also shows that the errors for features 1 and 5 are significantly worse 
than the errors for the other 3 features. The locations of these two features were 
often obscured by motion blur and specular reflections, hence the errors are due 
to inaccurate tracking and matching. 

The RMS errors for the « and y dimensions of features 2 to 4 range between 
1+0.62mm and 1.49+0.96mm. These are of the same order of magnitude as the 
noise we expect to be introduced by tracking and matching inaccuracies. This 
suggests that the CRPM’s estimate of what the compound cardio-respiratory 
motion would look like in the absence of noise may be adequate for it to be used 
for further analysis of that motion in these dimensions. 


Table 1 shows us that the spatial ranges of the middle 90% of the temporally- 
registered cardiac cycles of features 2 to 4 are still quite large (> 2.48mm) after 
subtracting the CRPM’s respiratory component (compare these ranges to the 
~20mm range of total motion shown in the left column of figure 2). These errors 
should be reduced before using the model to analyse the cardiac and respiratory 
motion components individually. One possible way to reduce them may be to 


Feature 2: Best-fitting model (n.=8, n =19, a=1.42) Registered cardiac cycles (n,=19, n,=8, 1=1.42) 
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Fig. 2. (Left column) These graphs show the fit of the CRPM to the zx, y and z 
trajectories of feature 2. (Right column) The red graphs show the spatial variations 
of feature 2’s cardiac cycles after subtracting the CRPM’s respiratory component and 
temporally registering them. The green graphs show the feature’s temporally registered 
cardiac cycles before the respiratory component was subtracted, for comparison. The 
thick central lines denote the median trajectories, and the error bars show 45% of 
the data either side of the median. The ratio of the mean height of the red error 
bars to the mean height of the green bars gives a measure of how much subtracting 
the CRPM’s respiratory component reduces the cardiac cycle variations. Minimising 
€CCR® will minimise this ratio and vice versa. The fact that this ratio is so high in 
the z dimension shows that the CRPM is fitting the z noise rather than accurately 
estimating the respiratory and cardiac motion. 


add a term to equation (2) for the CCR error, and to reformulate equation (3) 
to minimise this error in addition to the fitting error and respiratory curvature. 


4 Future Work 


It is infeasible for us to manually track enough features over multiple patients to 
be able to do a more comprehensive analysis at this stage. However, it is feasible 
for us to manually track a few features over a few cardiac cycles (one cycle is 
typically approximately 50 frames long). So the next step is for us to use short, 
manually-defined trajectories as prior information with which to parameterise 
the prior model of a rudimentary tracking algorithm. This should allow us to 
track salient features individually. 

Once we have sufficient data for these salient features, we will validate al- 
ternative cardio-respiratory motion decomposition methods. Accurate motion 
decomposition will allow us to develop better models of the motion of individual 
features. 

We will also evaluate models of the local, relative cardiac/respiratory motion 
of features. A good model of local, relative feature motion will allow us to track 
highly deforming regions and to improve tracking robustness by tracking salient 
features that are close to each other as a group. It will also provide prior infor- 
mation that we can use to infer the motion of features that lie in homogeneous 
regions near salient features. These new models will thus allow us to densely 
track patches on the myocardial surface. These patches can then be combined 
to form a reconstruction of the visible surface. 

To fully realise our myocardial reconstruction goals, the final task will be 
to solve the problem of automatically parameterising the prior models, so as to 
obviate the need to do any manual tracking for new videos. 


5 Conclusion 


We’ve shown how to estimate the cardiac and respiratory phases from intra- 
operative videos, and we’ve introduced a method for validating the results of 
any algorithm that decomposes cardio-respiratory motions. Furthermore, we’ve 
shown that the CRPM models the compound cardio-respiratory motion of our 
data set with adequate accuracy, but that its estimate of the respiratory mo- 
tion needs to be improved. A larger-scale study based on automatically-tracked 
features is needed before more general conclusions can be made. 
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